Degradation statistics

baRulho:quantifying habitat-induced degradation of (animal) acoustic signals

Published

October 2, 2024

Source code, data and annotation protocol found at https://github.com/maRce10/barulho_paper

1 Purposes

  • Estimate degradation metrics on re-recorded signals from playback experiment at Bosque de Tlalpan, Mexico City, 2019

  • Run statistical analyses

 

Load packages

Code
# install/ load packages
sketchy::load_packages(packages = c("remotes", "kableExtra", "knitr",
    "formatR", "rprojroot", "maRce10/warbleR", "ggplot2", "tidyr",
    "viridis", "corrplot", "brms", "ggdist", "cowplot", "cmdstanr",
    "maRce10/brmsish"))

source("~/Dropbox/R_package_testing/brmsish/R/extended_summary.R")
source("~/Dropbox/R_package_testing/brmsish/R/helpers.R")
Code
degrad_df <- read.csv("./data/processed/tlalpan_degradation_metrics_v01.csv")

degrad_params <- c("blur.ratio", "spectrum.blur.ratio", "envelope.correlation",
    "excess.attenuation", "signal.to.noise.ratio", "cross.correlation",
    "tail.to.signal.ratio", "tail.to.noise.ratio", "spectrum.correlation")

comp.cases <- complete.cases(degrad_df[, names(degrad_df) %in% degrad_params])

pca <- prcomp(degrad_df[comp.cases, names(degrad_df) %in% degrad_params],
    scale. = TRUE)

# add to data
degrad_df$pc1 <- NA
degrad_df$pc1[comp.cases] <- pca$x[, 1]
degrad_df$pc1.1m.rate <- degrad_df$distance
# plot rotation values by PC
pca_rot <- as.data.frame(pca$rotation[, 1:4])

pca_rot_stck <- stack(pca_rot)

pca_rot_stck$variable <- rownames(pca_rot)
pca_rot_stck$Sign <- ifelse(pca_rot_stck$values > 0, "Positive", "Negative")
pca_rot_stck$rotation <- abs(pca_rot_stck$values)

ggplot(pca_rot_stck, aes(x = variable, y = rotation, fill = Sign)) +
    geom_col() + coord_flip() + scale_fill_viridis_d(alpha = 0.7,
    begin = 0.3, end = 0.8) + facet_wrap(~ind)

1.1 Correlation among metrics

Raw metrics:

Code
cormat <- cor(degrad_df[, degrad_params], use = "pairwise.complete.obs")

rownames(cormat) <- colnames(cormat) <- degrad_params

cols_corr <- colorRampPalette(c("white", "white", viridis(4, direction = -1)))(10)

cp <- corrplot.mixed(cormat, tl.cex = 0.7, upper.col = cols_corr,
    lower.col = cols_corr, order = "hclust", lower = "number", upper = "ellipse",
    tl.col = "black")

Code
# sort parameters as in clusters for cross correlation
degrad_params <- degrad_params[match(rownames(cp$corr), degrad_params)]

1.2 Data description

  • 20 frequencies
  • 3 locations
  • 11520 test sounds
  • 160 sound treatment combinations
  • Sample sizes per location, transect and signal type
Code
agg <- aggregate(cbind(sound.id, treatment.replicates) ~ location +
    habitat.structure + distance, degrad_df, function(x) length(unique(x)))

agg$replicates <- agg$sound.id/agg$treatment.replicates

agg
location habitat.structure distance sound.id treatment.replicates replicates
1 closed 10 480 160 3
2 closed 10 480 160 3
3 closed 10 480 160 3
1 open 10 480 160 3
2 open 10 480 160 3
3 open 10 480 160 3
1 closed 30 480 160 3
2 closed 30 480 160 3
3 closed 30 480 160 3
1 open 30 480 160 3
2 open 30 480 160 3
3 open 30 480 160 3
1 closed 65 480 160 3
2 closed 65 480 160 3
3 closed 65 480 160 3
1 open 65 480 160 3
2 open 65 480 160 3
3 open 65 480 160 3
1 closed 100 480 160 3
2 closed 100 480 160 3
3 closed 100 480 160 3
1 open 100 480 160 3
2 open 100 480 160 3
3 open 100 480 160 3

1.3 Statistical analysis

Code
iter <- 10000
chains <- 4
priors <- c(prior(normal(0, 6), class = "b"), prior(cauchy(0, 4),
    class = "sd"))


# set frequency to mean-centered
degrad_df$frequency <- scale(degrad_df$frequency, scale = FALSE)

# set base level for factors
degrad_df$habitat.structure <- factor(degrad_df$habitat.structure,
    levels = c("open", "closed"))
degrad_df$frequency.modulation <- factor(degrad_df$frequency.modulation,
    levels = c("no_fm", "fm"))
degrad_df$amplitude.modulation <- factor(degrad_df$amplitude.modulation,
    levels = c("no_am", "am"))
degrad_df$duration <- factor(degrad_df$duration, levels = c("short",
    "long"))
degrad_df$location <- as.factor(degrad_df$location)
# degrad_df$ord.distance <- ordered(degrad_df$distance)

set.seed(123)

cmdstanr::set_cmdstan_path("~/Documentos/cmdstan/")

# to run within-chain parallelization
mod_pc1 <- brm(formula = pc1 ~ frequency * habitat.structure + frequency.modulation *
    habitat.structure + amplitude.modulation * habitat.structure +
    duration * habitat.structure + mo(distance) + (1 | location) +
    (1 | treatment.replicates), data = degrad_df, prior = priors,
    iter = iter, chains = chains, cores = chains, control = list(adapt_delta = 0.99,
        max_treedepth = 15), file = "./data/processed/regression_model_pc1.RDS",
    file_refit = "always", backend = "cmdstanr", threads = threading(10))

mod_blurratio <- brm(formula = blur.ratio ~ frequency * habitat.structure +
    frequency.modulation * habitat.structure + amplitude.modulation *
    habitat.structure + duration * habitat.structure + mo(distance) +
    (1 | location) + (1 | treatment.replicates), data = degrad_df,
    prior = priors, iter = iter, chains = chains, cores = chains,
    control = list(adapt_delta = 0.99, max_treedepth = 15), file = "./data/processed/regression_model_blur_ratio.RDS",
    file_refit = "always", backend = "cmdstanr", threads = threading(10))
Code
my.viridis <- function(...) viridis(alpha = 0.5, begin = 0.3, end = 0.7,
    ...)


effects <- c("habitat structure", "frequency modulation", "amplitude modulation",
    "frequency", "duration", "frequency:habitat structure", "habitat structure:duration",
    "habitat structure:amplitude modulation", "habitat structure:frequency modulation")
extended_summary(read.file = "./data/processed/regression_model_pc1.RDS",
    n.posterior = 1000, fill = "orange3", trace.palette = my.viridis,
    remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
        "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
        "durationlong"), gsub.replacement = c("", "habitat structure",
        "amplitude modulation", "frequency modulation", "duration"),
    )

1.4 regression_model_pc1

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 6) Intercept-student_t(3, -0.3, 2.5) sd-cauchy(0, 4) sigma-student_t(3, 0, 2.5) simo-dirichlet(1) pc1 ~ frequency * habitat.structure + frequency.modulation * habitat.structure + amplitude.modulation * habitat.structure + duration * habitat.structure + mo(distance) + (1 | location) + (1 | treatment.replicates) 10000 4 1 5000 42 (0.0021%) 0 1795.606 3839.76 401962879
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
frequency 0.157 0.127 0.187 1.001 2607.552 4448.647
habitat structure 1.086 0.990 1.181 1 10835.391 12469.973
frequency modulation 0.050 -0.129 0.227 1.002 1795.606 4018.020
amplitude modulation 0.268 0.090 0.442 1.002 2014.067 3839.760
duration -0.141 -0.320 0.033 1.002 2168.261 4800.424
frequency:habitat structure 0.047 0.031 0.063 1 19362.470 15343.730
habitat structure:frequency modulation -0.087 -0.180 0.005 1 15136.276 15435.642
habitat structure:amplitude modulation -0.007 -0.100 0.085 1 15877.418 14414.394
habitat structure:duration -0.003 -0.097 0.092 1 15285.902 13945.763

Code
extended_summary(read.file = "./data/processed/regression_model_blur_ratio.RDS",
    n.posterior = 1000, fill = "orange3", trace.palette = my.viridis,
    remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
        "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
        "durationlong"), gsub.replacement = c("", "habitat structure",
        "amplitude modulation", "frequency modulation", "duration"),
    effects = effects)

1.5 regression_model_blur_ratio

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 6) Intercept-student_t(3, 0.1, 2.5) sd-cauchy(0, 4) sigma-student_t(3, 0, 2.5) simo-dirichlet(1) blur.ratio ~ frequency * habitat.structure + frequency.modulation * habitat.structure + amplitude.modulation * habitat.structure + duration * habitat.structure + mo(distance) + (1 | location) + (1 | treatment.replicates) 10000 4 1 5000 100 (0.005%) 0 1905.433 4789.244 299027454
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
habitat structure 0.029 0.026 0.033 1 18137.687 14919.455
frequency modulation 0.065 0.058 0.072 1.002 2043.117 4835.342
amplitude modulation 0.023 0.016 0.030 1.002 1905.433 4789.244
frequency 0.001 0.000 0.002 1.001 2413.543 5026.225
duration 0.000 -0.006 0.008 1.003 2110.331 4830.679
frequency:habitat structure 0.003 0.002 0.003 1 20080.498 12256.359
habitat structure:duration 0.005 0.001 0.009 1 25094.631 15713.549
habitat structure:amplitude modulation 0.002 -0.001 0.006 1 25475.014 14808.890
habitat structure:frequency modulation -0.020 -0.024 -0.017 1 27043.019 14561.904

Code
extended_summary(read.file = "./data/processed/regression_model_excess_attenuation.RDS",
    n.posterior = 1000, fill = "orange3", trace.palette = my.viridis,
    remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
        "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
        "durationlong"), gsub.replacement = c("", "habitat structure",
        "amplitude modulation", "frequency modulation", "duration"),
    effects = effects)

1.6 regression_model_excess_attenuation

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 6) Intercept-student_t(3, 81.7, 31.2) sd-cauchy(0, 4) sigma-student_t(3, 0, 31.2) simo-dirichlet(1) excess.attenuation ~ frequency + frequency.modulation + amplitude.modulation + duration + habitat.structure + mo(distance_f) + (1 | location) + (1 | treatment.replicates) 20000 4 1 10000 5 (0.00012%) 0 3308.454 5800.994 1741877593
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
habitat structure 7.486 7.072 7.898 1 42076.716 29423.336
frequency modulation -3.427 -4.985 -1.869 1.002 3308.454 7139.187
amplitude modulation -3.295 -4.841 -1.686 1.001 3470.592 5800.994
frequency 13.703 12.116 15.299 1.001 3473.790 7645.888
duration 0.792 -0.824 2.390 1.001 3656.833 6971.910
NA NA NA NA NA NA NA
NA.1 NA NA NA NA NA NA
NA.2 NA NA NA NA NA NA
NA.3 NA NA NA NA NA NA

Code
# extended_summary(read.file =
# './data/processed/regression_model_signal_to_noise_ratio.RDS',
# n.posterior = 1000, fill = 'orange3', trace.palette =
# my.viridis, remove.intercepts = TRUE, highlight = TRUE,
# gsub.pattern = c('b_', 'habitat.structureclosed',
# 'amplitude.modulationam', 'frequency.modulationfm',
# 'durationlong'), gsub.replacement = c('', 'habitat structure',
# 'amplitude modulation', 'frequency modulation', 'duration'),
# effects = effects)

extended_summary(read.file = "./data/processed/regression_model_tail_to_signal_ratio.RDS",
    n.posterior = 1000, fill = "orange3", trace.palette = my.viridis,
    remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
        "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
        "durationlong"), gsub.replacement = c("", "habitat structure",
        "amplitude modulation", "frequency modulation", "duration"),
    effects = effects)

1.7 regression_model_tail_to_signal_ratio

priors formula iterations chains thinning warmup diverg_transitions rhats > 1.05 min_bulk_ESS min_tail_ESS seed
1 b-normal(0, 6) Intercept-student_t(3, -6, 8.6) sd-cauchy(0, 4) sigma-student_t(3, 0, 8.6) simo-dirichlet(1) tail.to.signal.ratio ~ frequency + frequency.modulation + amplitude.modulation + duration + habitat.structure + mo(distance_f) + (1 | location) + (1 | treatment.replicates) 20000 4 1 10000 26 (0.00065%) 0 3803.126 7149.158 1695545447
Estimate l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
habitat structure 3.270 3.099 3.439 1 44650.292 29637.846
frequency modulation -0.692 -1.297 -0.085 1.001 3901.160 7710.797
amplitude modulation 0.659 0.049 1.257 1 3803.126 8310.054
frequency 1.622 1.025 2.211 1 4186.279 8273.165
duration -0.400 -1.006 0.202 1.002 3961.338 7149.158
NA NA NA NA NA NA NA
NA.1 NA NA NA NA NA NA
NA.2 NA NA NA NA NA NA
NA.3 NA NA NA NA NA NA

 

2 Takeaways

 


 

Session information

R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.4 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       

time zone: America/Costa_Rica
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] brmsish_1.0.0       cmdstanr_0.8.1.9000 cowplot_1.1.3      
 [4] ggdist_3.3.2        brms_2.21.0         Rcpp_1.0.13        
 [7] corrplot_0.92       viridis_0.6.5       viridisLite_0.4.2  
[10] tidyr_1.3.1         ggplot2_3.5.1       warbleR_1.1.32     
[13] NatureSounds_1.0.4  seewave_2.2.3       tuneR_1.4.7        
[16] rprojroot_2.0.4     formatR_1.14        knitr_1.48         
[19] kableExtra_1.4.0    remotes_2.5.0      

loaded via a namespace (and not attached):
  [1] bitops_1.0-8         pbapply_1.7-2        gridExtra_2.3       
  [4] inline_0.3.19        testthat_3.2.1.1     sandwich_3.1-0      
  [7] rlang_1.1.4          magrittr_2.0.3       multcomp_1.4-25     
 [10] matrixStats_1.3.0    compiler_4.4.1       loo_2.8.0           
 [13] reshape2_1.4.4       systemfonts_1.1.0    vctrs_0.6.5         
 [16] stringr_1.5.1        arrayhelpers_1.1-0   pkgconfig_2.0.3     
 [19] crayon_1.5.3         fastmap_1.2.0        backports_1.5.0     
 [22] labeling_0.4.3       fontawesome_0.5.2    utf8_1.2.4          
 [25] rmarkdown_2.28       ps_1.7.7             purrr_1.0.2         
 [28] xfun_0.47            jsonlite_1.8.8       highr_0.11          
 [31] tidybayes_3.0.6      uuid_1.2-1           parallel_4.4.1      
 [34] R6_2.5.1             stringi_1.8.4        StanHeaders_2.32.10 
 [37] brio_1.1.5           estimability_1.5.1   rstan_2.32.6        
 [40] zoo_1.8-12           bayesplot_1.11.1     Matrix_1.7-0        
 [43] splines_4.4.1        tidyselect_1.2.1     rstudioapi_0.16.0   
 [46] abind_1.4-5          yaml_2.3.10          dtw_1.23-1          
 [49] codetools_0.2-20     processx_3.8.4       curl_5.2.2          
 [52] pkgbuild_1.4.4       plyr_1.8.9           lattice_0.22-6      
 [55] tibble_3.2.1         withr_3.0.1          bridgesampling_1.1-2
 [58] posterior_1.6.0      coda_0.19-4.1        evaluate_0.24.0     
 [61] signal_1.8-1         survival_3.7-0       sketchy_1.0.3       
 [64] proxy_0.4-27         RcppParallel_5.1.8   xml2_1.3.6          
 [67] pillar_1.9.0         tensorA_0.36.2.1     packrat_0.9.2       
 [70] checkmate_2.3.2      stats4_4.4.1         distributional_0.4.0
 [73] generics_0.1.3       RCurl_1.98-1.16      rstantools_2.4.0    
 [76] munsell_0.5.1        scales_1.3.0         xtable_1.8-4        
 [79] glue_1.7.0           emmeans_1.10.3       tools_4.4.1         
 [82] xaringanExtra_0.8.0  mvtnorm_1.2-5        grid_4.4.1          
 [85] ape_5.8              QuickJSR_1.3.1       colorspace_2.1-1    
 [88] nlme_3.1-165         cli_3.6.3            svUnit_1.0.6        
 [91] fansi_1.0.6          svglite_2.1.3        Brobdingnag_1.2-9   
 [94] dplyr_1.1.4          V8_4.4.2             gtable_0.3.5        
 [97] fftw_1.0-8           digest_0.6.37        TH.data_1.1-2       
[100] farver_2.1.2         rjson_0.2.22         htmlwidgets_1.6.4   
[103] htmltools_0.5.8.1    lifecycle_1.0.4      MASS_7.3-61